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The measurement of photon-number statistics of fields composed of photon pairs, generated in 
spontaneous parametric down-conversion and detected by an intensified CCD camera is described. 
Final quantum detection efficiencies, electronic noises, finite numbers of detector pixels, transverse 
intensity spatial profiles of the detected beams as well as losses of single photons from a pair are taken 
into account in a developed general theory of photon-number detection. The measured data provided 
by an iCCD camera with single-photon detection sensitivity are analyzed along the developed theory. 
Joint signal-idler photon-number distributions are recovered using the reconstruction method based 
on the principle of maximum likelihood. The range of applicability of the method is discussed. The 
reconstructed joint signal-idler photon-number distribution is compared with that obtained by a 
method that uses superposition of signal and noise and minimizes photoelectron entropy. Statistics 
of the reconstructed fields are identified to be multi-mode Gaussian. Elements of the measured as 
well as the reconstructed joint signal-idler photon-number distributions violate classical inequalities. 
Sub-shot-noise correlations in the difference of the signal and idler photon numbers as well as partial 
suppression of odd elements in the distribution of the sum of signal and idler photon numbers are 
observed. 

PACS numbers: 42.65.Lm,42.50.Ar 
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I. INTRODUCTION 

Light generated in the process of spontaneous paramet- 
ric down-conversion (SPDC) is emitted in photon pairs 
Photons comprising one photon pair are strongly 
quantum correlated (entangled). Entanglement of pho- 
tons in a pair has been used in many experiments that 
have provided a deep insight into the laws of quantum 
mechanics [3, 's^. Among others, the measured violation 
of Bell's inequalities ruled out neoclassical local hidden- 
variables theories. Photon pairs have also found their 
way topractical applications, e.g., in quantum cryptog- 
raphy [j], measurement of ultrashort time delays, or ab- 
solute measurements of detection quantum efficiencies 
0, . These experiments utilize photon fields that con- 
tain only one photon pair in a measured time window 
with a high probability. 

There have been experiments (teleportation, measure- 
ment of GHZ correlations, etc.) measuring triple and 
quadruple coincidence counts caused by fields containing 
two photon pairs in a time window given by an ultrashort 
pump pulse. However, states used for such experiments 
contain a very low fraction of states with two photon 
pairs in comparison with the fraction belonging to the 
state with one photon pair and the vacuum state. The 
reason is to eliminate the influence of three and more- 
than-three photon-pair states to the considered experi- 
mental setups. Measurements done in such setups have 
to be conditional and they require long data-acquisition 



times. 

The use of more powerful pump pulses as well as devel- 
opment of materials with higher values of x*-^' suscepti- 
bilities have opened the way to generate fields containing 
many photon pairs originating in one pump pulse. For 
such fields, a joint signal-idler photon-number distribu- 
tion is the main characteristic that determines the ex- 
perimental results utilizing these fields. Determination 
of photon-pair statistics is important also for weak cw 
fields provided that they are detected in long-time de- 
tection windows [7]. In this case photon-pair statistics 
have been identified to be Poissonian after eliminating 
dead-time detection effects Q- 

Returning back to more intense fields, recent exper- 
iments f8!-fl5!| (and references therein) are even able to 
provide experimental joint signal-idler photoelectron dis- 
tributions of twin beams containing up to several thou- 
sands of photon pairs. As for detectors, weaker fields 
containing up to ten photons can be measured by special 
single-photon avalanche detectors (VLPC) [l^, hybrid 
photo- multipliers [TtI . [isj . super-conducting bolometers 
II9I or time-multiplexed fiber-optics detection loops f20l - 
l25j |. Intensified CCD cameras [ll|, [2^ [23| can in prin- 
ciple capture states with hundreds of photons. Ultra- 
sensitive photodiodes with their linear response and very 
low level of noise are suitable for the detection of states 
with hundreds or even better thousands of photons. A 
special method utilizing precisely attenuated beams has 
also been suggested and developed [1^, [2§| . It allows to 
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resolve photon numbers even in the measurement based 
on single-photon sensitive avalanche photodiodes. We 
note that also the well-known homodyne detection has 
been found useful in the determination of intensity cor- 
relations of twin beams [10, IMj . 

All these approaches give experimental photoelectron 
distributions obtained by detectors with finite quantum 
detection efficiencies. While silicon PIN photodiodes or 
back-illuminated CCD cameras can offer detection effi- 
ciencies close to unity, their internal noise prevents their 
use in the single-photon regime. On the other hand, de- 
tectors with large internal gain, like iCCD cameras, EM- 
CCD cameras or avalanche photo-detectors allow single- 
photon sensitivity or even photon-number resolution by 
effectively decreasing their noise. A price for this sen- 
sitivity is paid, however, in the form of lower quantum 
efficiencies. ICCD cameras are a good trade-off in this re- 
spect. Their level of noise is low, but not negligible. On 
the other hand, quantum detection efficiencies around 
20 % are sufficient enough to profit from their low level 
of noise. 

Once we know quantum detection efhciency and the 
level of noise, we can reconstruct the field in front of a 
detector. The usual and physically-motivated approach 
is based on the assumption of the character of the recon- 
structed field. Working with photon pairs, we can nat- 
urally assume that the reconstructed field is composed 
of certain number of independent modes containing pho- 
ton pairs and small additional noise in the form of sin- 
gle photons. Using this picture, a multi-mode theory of 
signal and noise tailored specially for paired fields can 
be applied (see [s^, for the spontaneous process and 
(34 - . 36] for the stimulated process). It can be accompa- 
nied with the principle of minimum entropy to get the 
reconstructed field. As an alternative, one may rely on a 
mathematically based method that uses the maximum- 
likelihood principle. In the framework of this method, the 
reconstructed field is reached as a steady point accessible 
by an iteration procedure. As a final step in characteriza- 
tion of the fields, joint signal-idler quasi-distributions of 
integrated intensities may be reached using the recon- 
structed joint signal-idler photon-number distributions 

[Sill. 

Here, we pay attention to the determination of a joint 
signal-idler photon-number distribution beyond a non- 
linear crystal using an iCCD camera as a tool resolving 
photon numbers. The method of maximum likelihood is 
applied. It allows to deal with even more difficult experi- 
mental conditions like those reached when more than one 
photon can be registered in a single pixel. Transforma- 
tion matrices describing details of the detection process 
and being an important ingredient of the iteration pro- 
cedure are derived under several conditions. The recon- 
structed fields are compared with those obtained by the 
method of superposition of signal and noise. 

The paper is organized as follows. Sec. II contains a 
general model describing a photon-number-resolving de- 
tection device. Sec. Ill is devoted to the description of 




FIG. 1; Scheme of the considered model. Photon pairs are 
generated in a nonlinear crystal NLC. Virtual beam-splitters 
BSs and BSj describe possible losses of one or both photons 
from a pair before they are detected. Signal (idler) photons 
propagate through a 1 x Ns (1 x Nj) multi-port and are de- 
tected by single-photon detectors Dsi , Ds2, • • • , '-'Sns i^ii , 
D/2, . . . , D/„^). Signals from the detectors are registered in 
a coincidence-detection device. 

photon-number-resolving detection by an iCCD camera 
under real experimental conditions. The role of inhomo- 
geneous transverse profiles of the detected fields is dis- 
cussed in Sec. IV. In Sec. V, the iteration procedure of 
the maximum-likelihood method is explained and used 
to recover joint signal-idler photon-number distributions. 
Subsec. VA is devoted to nonclassical characteristics of 
the emitted fields. Statistics of the fields are discussed in 
Subsec. VB in which also the problem of reconstruction of 
more intense fields is addressed. Comparison of the re- 
constructed fields obtained by the maximum-likelihood 
method and the method of superposition of signal and 
noise is provided in Sec. VI. Sec. VII brings conclusions. 
The formula for an effective quantum detection efficiency 
is derived in Appendix A. 

II. PROBABILITIES OF 
MULTI-PHOTON-COINCIDENCE COUNTS IN 
AN ARRAY OF SINGLE-PHOTON DETECTORS 

The measurement of joint signal-idler photon-number 
distribution can be in general described using the scheme 
shown in Fig. [T] |38l - l4l| . Photon pairs occurring in the 
output plane of a nonlinear crystal NLC propagate to- 
wards photon-number-resolving detection devices placed 
in the paths of the signal and idler fields. One or both 
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photons from a pair may be lost before they reach their 
detection devices. Geometric filtering (one photon from 
a pair is not steered to the detector area), reflections on 
optical elements in the experimental setup or absorption 
of a photon along its path to a detection device repre- 
sent possible reasons. We describe this effect by two 
beam-splitters BSg and BS/ [l^l placed in the signal- 
and idler-field paths, respectively. We model a photon- 
number-resolving detection device as a multi-port 1 x iV 
[43} followed by N single-photon detectors. This descrip- 
tion holds also in the special case when an intensified 
CCD camera is used [2g|. We note that detectors able 
to resolve directly photon numbers to some extent have 
been constructed [ll, From practical viewpoint, 

also detectors using time multiplexing (reached, e.g., in 
fiber optics) and one or two single-photon detectors are 
promising (21I [22l . [25| . Photon- number-resolving detec- 
tion in all these devises can be described in the framework 
of the presented general model. 

We assume that the signal and idler fields in the out- 
put plane of the nonlinear crystal NLC are described by 
the following statistical operator psi written in the Fock 
basis: 
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PSI ^ ^ ^ p{ns,ni)\ns)ss{ns\<S)\ni)ii{ni\; (1) 

ns—O ni—Q 

the symbol p{ns,ni) stands for the joint signal-idler 
photon-number distribution. 

Statistical operator p^j appropriate for the signal and 
idler fields in front of the detection devices can be ex- 
pressed as t3&] : 
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The symbols Rs and Ri {Ts and Ti) denote intensity re- 
flectivities (transmissivities) of the beam-splitters in the 
corresponding path. 

We assume a multi-port f x Ns {IxNj) followed by Ns 
(Nj) single-photon detectors with quantum efficiencies 
VSj ivij) s-nd dark-count rates ds^ i^ij) the signal 
(idler) path. Detection of a photon in the fc-th detector 
is described by the following detection operator Di^ [3^ : 
00 

A. = J2^[l-{l-V^X]+d^A^~V^X}\n)kk{n\, 



i = S,I. 



(3) 



On the other hand, detection operator Df° corresponds 
to the case when no detection has occurred: 
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(4) 



The effect of 'splitting' photons in the signal field in 
a 1 X Ns multi-port can be described by the relation 
as = ^2^=1 ^Sj aSj , where as is the annihilation opera- 
tor in the signal field entering the multi-port, whereas 
the annihilation operator ds describes a field at the i-th 
multi-port output. Symbol ts^ stands for an amplitude 
transmissivity of a photon from the input to the j-th 
output. The 1 X Nj multi-port in the idler- field path is 
described similarly and the symbol tj. then refers to an 
amplitude transmissivity of a photon from the input to 
the j-th multi-port output. 

The probability Csd jd that given eg detectors in the 
signal field and given c/ detectors in the idler field detect 
a photon whereas the rest of detectors does not register 
a photon is determined as the quantum mean value: 
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The symbol S (/) denotes the set of all signal-field (idler- 
field) detectors S = {Si, . . . , Sns} {I = {h, ■ ■ ■ ,Ini})- 
The set {I^) contains signal-field (idler-field) detec- 
tors that have registered a photon. Symbol Tr stands for 
an operator trace. 

Using the statistical operator pgj given in Eq. ([2]) the 
probability C5D of a multi-coincidence count defined 
in Eq. ^ is obtained in the form: 
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We now consider two symmetric multi-ports (tg^ = 
ts2 = ■ ■ ■ tstis — '^s = ^/V^s, tj^ = ti^ — . . . — = 
tl — 1/ y/TT]) and detectors endowed with the same char- 
acteristics in the signal and idler fields (775^ — rjg^ — 
• • • = "HSns = ^s, ^ ds^ = ■ ■ ■ = ds^^ = ds, Vh = 
VI2 = • • ■ = VInj = Vi, dj^ ^ di^ = . . . = dj^^ = dj). 
Then the probability f'^^''^'{cs,ci) of having eg detec- 
tions somewhere at Ns signal detectors and c/ detections 
somewhere at Nj idler detectors can be expressed as: 

r-"^(cs,c,) = (^;)(^^)c,.^,.. (7) 

Using the expression for Cgo jd provided in Eq. ([6]) we 
arrive at the relation: 
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Ti (ji — TiTji) determines the probability that a photon 
is registered at some of the detectors. 

If the number of photons detected by the camera is 
much lower than the number of pixels detecting the 
overall field with a non-negligible probability, the lim- 
its Ns — > 00 and Nj — > 00 are appropriate. When 
determining these limits, the overall noise levels Ds and 
Dj are kept constant {Ds = Nsds, Dj = Njdi). The 
coefficients K defined in Eq. ^ then considerably sim- 
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We note that the following relations have been used when 
deriving Eq. ([T0| : 
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iV > 0;0° = 1; 
symbol denotes positive integer numbers. 



III. PHOTON-NUMBER DETECTION UNDER 
REAL EXPERIMENTAL CONDITIONS 

In our typical experiment (see Fig. [5]) we define three 
regions-of-interest on the iCCD detection photocathode: 
one for collecting signal photons, one for counting idler 
photons and the third one that serves for monitoring of 
the dark noise in the experiment. To achieve higher data 
collection rates we use hardware binning of several pixels 
to a single macro-pixel. The signal and idler regions typ- 
ically contain about one thousand of macro-pixels that 
give an information about photons detection. This means 
that a finite number of (macro-)pixels may play an im- 
portant role depending on intensity of the impinging field 
and the general form of the transfer matrix K^'^' {ci,ni) 
in Eq. ^ should be preferably used. However, evalua- 
tion of a transfer matrix K for larger numbers of photons, 
photoelectrons (registered photons) and (macro-)pixels is 
difficult because a high extended precision in the evalua- 
tion of the sum occurring in Eq. 1^ is required. For in- 
stance, if fields having up to 1000 photons are measured, 
from 2 to 3 hundred significant decimal digits are needed 
in the evaluation of the sums in Eq. ^ under conditions 
considered below. This is time demanding and that is 
why we present several alternative ways how to handle 
the problem under specific conditions. 

First, we rewrite the relation between the measured 
frequencies f^s,Ni ^^le photon-number distribution 
p in a general form: 
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where the general transformation matrices G^'^' {ci,ni) 
for i = S, I have been introduced. 
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In a real experimental setup, there are non-negligible 
losses (described by intensity transmissivities T5 and Tj) 
before a field arrives to the photocathode of the iCCD 
camera. As a result an average number of photons in 
the input to the camera is lower compared to the aver- 
age number of photons in the output plane of the crys- 
tal. This may make a numerical evaluation of the matrix 
K"^'^' in Eq. © faster due to lower dimensions of this 
matrix. In this case the transfer matrix G*'^* (q, rii) can 
be rewritten as a product of two matrices: 

Hi 

G^'^'(c„nO = i^*''^'(Q,m)i^o(m,n,). (12) 

The matrix Kq introduced in Eq. ([T^ describes the 
Bernoulli distribution with transmissivity Tf. 



Ko{m,ni) = 



(13) 



Evaluation of the matrix K'^^^^ using Eq. ([9]) is then done 
assuming r,; = 77^. 

If numbers of photons in a detected field are too high 
preventing from the application of Eq. ^ (technical rea- 
sons in the evaluation) we can proceed as follows. The 
measured field first undergoes losses described by the in- 
tensity transmissivity Ti before impinging on the cam- 
era. In the next step each photon present in the region- 
of-interest of the camera 'registers' itself in one (macro- 
)pixel. The probability 7*'^* (77125 ™i) that mi photons 
register in m2 (macro-)pixels assuming the overall num- 
ber of (macro-) pixels to be Ni is given by permutations 
with repetition: 



7*'^'(m2,mi) 
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N, +7711-1 

N,, - 1 



m2 < Ni. (14) 



In this case, m2 (macro-)pixels are exposed by the field 
and the probability of detections (c; < 7712) is given 
by the matrix K'^'°°{ci,m2) written in Eq. (flUl) . i.e. as 
if there is an infinite number of (macro-) pixels in the 
camera. The matrix G^'^^ {ci,ni) then takes its final ap- 
proximative form: 

min(mi,A''i) m 

G^'^'(c„ri,) = E E K''^ic^,m2) 

m2— rni— 

X 7*'^'(7772,mi)i^o(mi,77,); (15) 



the matrix Kq is defined in Eq. (|13p . 

It has been assumed in the derivation of Eq. ([TSl) that 
each of the exposed 777,2 (macro)pixels contains only one 
photon (see the limit N — > 00). This approximation 
can be improved. If 77ii photons is registered at 7772 
(macro-)pixels, an average photon number occurring in 
one (macro-)pixel is 7711/7772. The average photon num- 
ber 7771/7772 greater than one leads to a higher probability 



of detection. This increase of detection probability can 
be modelled by an effective increase of detection quan- 
tum efiiciency (see Appendix A). The improved matrix 
G^'^' {ci,ni) can then be determined along the relation: 



G*'^'(c„ 77,) = r*''^'(c„7n)Ko(m,7 



(16) 
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On the other hand weak detected fields allow the fol- 
lowing simplification. If the maximum number Ci of 
counts is much less than the number of (macro-) pixels 
Ni the expression for matrix K'^'^^ in Eq. ^ can be suc- 
cessfully approximated using the relation (1 -I- x)"' « 
exp(a;77i) for \x\ <C 1. We then arrive at the matrix 
G^''^' (ci,ni) in the form: 



G''^'(c„ 77,) = g i^:,^' (c„ 777)^0 (m,".) 



(18) 
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and 



Kj^^{c,,m) = (^^^)(l-d.)^'-^'(l-7y,)"-^- 



di{l - tli) + Vi 



N 



(19) 



The expression in Eq. has a simple interpretation: 
m — Ci impinging photons are not registered with proba- 
bility 1 — r]i per photon. There occur Ci counts given ei- 
ther by impinging photons with probability mrji/Ni per 
photon or by dark counts with probability di{l — rji). 
Ni — Ci (macro-)pixels cannot feel dark counts with prob- 
ability 1 — di per (macro-) pixel. 



IV. INHOMOGENEOUS TRANSVERSE 
INTENSITY PROFILE OF A DETECTED BEAM 

If the intensity transverse profile of a beam imping- 
ing on an iCCD camera is inhomogeneous, we can di- 
vide (macro-) pixels of the camera into Mi groups assum- 
ing the same level of illumination of (macro-)pixels be- 
longing to one group. A j-th group of (macro-)pixels is 
characterized by probability Ti . that a photon present in 
beam 7(7 = 5*, /) impinges on one (macro-)pixel from this 
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group, number Vi. of (macro-) pixels, quantum detection 
efRciency rji- dark-count rate di - , and mean number pL^. 
of photons impinging on one (macro-)pixel. It holds that 

'^ij^ij = Vi and I]jJi = The probability 
that a photon reaches one (macro-)pixel of the j-th group 
is linearly proportional to the mean number [ii. of pho- 
tons coming to this (macro-)pixel and can be expressed 
as: 



■■iV,; 



(20) 



The average mean photon number /if™"^ is given as 

A transformation matrix /iT''^' (cj, n^) that generalizes 
the matrix /sT*'^* occurring in Eq. ^ and determines the 
probability of counts caused by rii photons coming to 
the camera is given by the following Mj-dimensional con- 
volution of matrices K^-^^j written in Eq. ^ and charac- 
terizing the detection in the j-th group of (macro-) pixels: 



if''^'(c„n,) = 



Mi rii 
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X ni\ J]^ 



(21) 



The matrix K^-^^ occurring in Eq. (|2T|) can be rewrit- 
ten into the following form if the matrices K^'^^-i are ex- 
pressed using the relation in Eq. ([9|): 
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(22) 



If the number of (macro-)pixels is sufficiently large 
compared to the number of impinging photons, consider- 
ation of the following limit is useful. In this limit — > oo 
for j = 1, . . . , Mi assuming Vi Ti- [probability that a pho- 
ton is detected in the j-th group of (macro-)pixels] to be 
constant. Also di Vi- = Di [overall dark-count rate of 



all (macro-) pixels in the j-th group] is assumed to be 
constant. Then the expression in Eq. ([2^ simplifies and 
leaves the matrix fC'^' in the form: 
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The expression in Eq. (j23p has a simple interpretation: Ui 
photons impinging on the camera generates Ij counts in 
a j-th group of (macro-)pixels with probability Ti.Vi.rji. 
per photon and (cj — Ij) counts come from dark counts 
occurring in the j-th group of (macro-)pixels. The re- 
maining Hi — X]fc=i ^fc photons are not registered with 
probability 1 - Y.jtiinj'^zjVij) per photon. 

Provided that the maximum number of counts Cj in a 
j-th group of (macro-)pixels is much less than the num- 
ber of (macro-)pixels Vi^ in this group the approximate 
relation (1 -I- x)"' w exp(nia;) for \x\ 1 enables to re- 
arrange the formula in Eq. as follows: 
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(24) 



The expression in Eq. ()24p can be interpreted similarly 
as the formula in Eq. (jl9l) . If Ci counts occur after rii 
photons enter the camera, Ui — Ci photons is not regis- 
tered with probability 1 — X^fcf^i '''ik'^ikVik P^r photon. In 
a j-th group of (macro)pixels Vi . — Cj (macro-)pixels do 
not count a photon with probability 1 — c?i . per (macro- 
)pixel (dark counts have to be 'eliminated'). Finally Cj 
(macro-)pixels detect a photon either due to a successful 
registration of a photon taken from Ui incident photons 
with probability Ti.iji. per photon or owing to a dark 

count with probability c?i .(l — ^ifc^u ) dark 

count occurs if there is no detection caused by an im- 
pinging photon). 
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If the number of counts ci registered by an iCCD cam- 
era is low and the number of groups of (macro-)pixels 
is greater, a useful alternative expression for the trans- 
fer matrix JsT*'^' (q, n^) given in Eq. (PT|) can be derived 
directly from Eq. ([6]) (T^ = 1 and Ri = are assumed, 



Mi a 

HE 

X -^4,{ci,...xjv,j(Q,rij), i^S,I, (25) 

where the coefficient iiri_{ci,...,CM.}(ci, fii) gives the prob- 
ability that Cj counts have occurred in the j-th group 

of (macro-)pixels; J2jti — Ci- It can be expressed as 
follows: 
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[Bj + Tiai{ki)rii<T,iki) + Tiai{k2)Vi<Ti{k2)]""' 

(1 - diai{ki)){^ - dia,{k2)) 
Ci-2 Ci-1 Ci 

(-ir-^^E E E 

fci— 1 A;2 — /ci+l A;3— ^2 + 1 

3 1 
0» + Y.m=li'^i'^iikm)'n^<T,{k^)) 



nm=l(l"C?»<T.(fc„)) 
[0i + Em=l(TirT,(fc„)»?icT.(fc„))]"* 



(26) 



The symbol <di introduced in Eq. (|26|) denotes the prob- 
ability that a photon is not registered by the camera; 
i.e. 

Mi 

e, = l-Y,ir^,'^^,V^,), l = S,I. (27) 
J = l 

The vector ai in Eq. (j26p is composed of elements 
(i = 5*, /) ; its j-th element gives the number of group 
of (macro-) pixels that registered the j-th click (j — 
1, . . . ,Ci). Thus, the first ci elements equal 1, the next 
C2 elements equal 2 and so on. 



V. RECONSTRUCTION OF THE JOINT 
SIGNAL-IDLER PHOTON-NUMBER 
DISTRIBUTION 

The probabilities (frequencies) /(c5,c/) are measured 
in the experiment and the relation in Eq. (jlip has to be 



inverted in order to obtain the joint signal-idler photon- 
number distribution pins, nj) beyond the nonlinear crys- 
tal. The relation in Eq. ([Tl]) together with the coeffi- 
cients G''^'(ci,ni) defined in Eq. ^ can be inverted 
under special conditions analytically. For instance, if 
Ds — Dj = the inversion relation is found using 
the 'convolution' of function / with the Bernoulli dis- 
tributions with efficiencies l/iTsrjs) and l/{Tjrij) that 
are greater than one. For more general method 
of direct matrix inversion has been elaborated [1^, . 
However, analytical approaches are not suitable for pro- 
cessing real experimental data (47j because of numerical 
instabilities and the occurrence of artifacts. On the other 
hand reconstruction algorithms have occurred to be suit- 
able for this task [21] . Such algorithms are able to find a 
reconstructed joint signal-idler photon-number distribu- 
tion Picc('T^5, JT-/) that matches the measured frequencies 
/(csjC/) in the best way with respect to a given crite- 
rion. Here, we consider the KuUback-Leibler divergence 
as a measure of the distance between the experimental 
data and data provided by the developed theory. The 
reconstructed joint signal-idler photon-number distribu- 
tion PrecinsjUi) minimizing the KuUback-Leibler diver- 
gence can then be found as a steady point of an iteration 
algorithm [Mlii]: 

p("+i)(n5,nj) = p^^\ns,nj) 

(28) 

The symbol p^'^\ns,ni) denotes a joint signal-idler 
photon-number distribution after an n-th step of the iter- 
ation, p^'^\ng,nj) is an arbitrary initial photon-number 
distribution. We note that each element of the initial 
photon-number distribution has to be nonzero in order 
to be considered in the iteration process. 

Convergence of the iteration process can be monitored 
using parameter S that gives the declination of the recon- 
structed photon-number distribution from the measured 
frequencies /(cs, c/) and is determined along the formula: 



n) 



E /('=s,c/) 



.cs,c/=0 



- G^'^^(CS,JS)G^'^^(C/,J7)P(")(J5,J7) 



1/2 



(29) 



Alternatively, also covariance C of the signal and idler 
photon numbers ng and ni derived for the joint photon- 
number distribution p, 



(AnsAn/) 



V((Ans)2)((An,)2)' 
An,; = Tii - (ni), i = 



(30) 
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FIG. 2: Scheme of the setup for detection of photon pairs. 
Fields composed of typically tens or hundreds of photon pairs 
are generated in a nonlinear crystal. Idler photons are re- 
flected on a mirror. Both signal and idler photons propagate 
through a frequency filter and are detected in an iCCD cam- 
era. 

oo oo 
ns—0 nj—0 

fc,/ = 0, 1,..., 

can be used as a useful indicator. The reason is that the 
initial photon-number distribution p^^^ is usually consid- 
ered without any correlation and the iteration process 
gradually reveals photon-number correlations present in 
the joint signal and idler field. We note that we have 
checked by numerical simulations that the reconstruc- 
tion algorithm cannot reveal correlations provided that 
the measured frequencies f(cs,ci) describe two indepen- 
dent fields. 

Here, we analyze three photon-number distributions 
obtained under different conditions using the experimen- 
tal setup shown in Fig. [2j Photon pairs are generated 
in a 5-mm long BBO crystal cut for a type I process 
{9 — 50 deg, (/) = 90 deg) pumped by ultrashort pulses de- 
livered by a cavity-dumped titanium-sapphire femtosec- 
ond laser at the wavelength of 840 nm followed by a 
third-harmonic generator. The pulses at the fundamental 
wavelength are about 150 fs long. The laser system runs 
at the repetition rate of 50 kHz and, after converting the 
840-nm beam to its third harmonic, it typically delivers 
pulses with the energy up to 45 nJ. Degenerate signal and 
idler photons occur at the cone layer behind the crystal 
and leave the crystal at the outer output half- angle of 13 
degrees. Photons in the idler field are reflected from a 
high-reflectivity mirror (^99 % at 560 nm) and impinge 
on an intensified CCD camera (Andor iStar 734). The 
camera has one megapixel resolution with 13 x 13 /im^ 
pixels but we use 8x8 hardware binning to increase the 
readout rate. The detection events are processed using 
our own software employing thresholding, event centroid 
finding and photon counting. The software is optimized 
to achieve maximum detection efficiency which we eval- 
uate to 23% near the degenerate wavelength of 560 nm. 
Three regions of interest are defined in the field of view 
of the camera: the first one is for the signal field, the 
second one for the idler field, and the last region serves 
for reference measurements of the noise level. The whole 
field of the camera is filtered by a high-transmittance 



high-pass filter blocking light below 490 nm and an in- 
terference filter of 14 nm FWHM centered at 560 nm. 
The interference filter selects nearly frequency degener- 
ate photon pairs. Since a single run of data acquisition 
usually takes several hours, the laser intensity is actively 
stabilized (the intensity noise lays below 0.3 % rms) using 
feedback loop and polarization attenuator. The intensi- 
fler of the camera is synchronously gated (gate duration 
equals 5 ns) by cavity-dumper trigger pulses to minimize 
the noise from the laboratory. 

In the experimental setup, histograms f{cs, c/) of pho- 
toelectron numbers have been taken under two different 
intensity conditions. In the first case (a), the measure- 
ment has been performed for lowest signal and idler inten- 
sities allowed by the setup. The limiting intensities are 
given by noise of the camera and stray light from the lab- 
oratory. The second case (b) represents a typical result 
obtained under most suitable conditions. The third case 
(c) corresponds to the measurement done with greater 
signal and idler intensities and the histogram /(cs,c/) 
has been obtained by summing up five neighbor frames 
together. We thus have three representative data sets 
with mean photoelectron numbers equal to 1.2, 8.6, and 
43. The corresponding histograms f{cs, cj) are shown in 
Fig. [21 Covariances of the photoelectron numbers cs and 
c/ described by the histograms /(eg, c/) plotted in Fig. [3] 
are in turn 23.8, 21.4, and 21.1 %. This corresponds to 
the expected overall detection efficiencies Tsrjs and Tjrji 
around 20 % and the low level of single-photon noise. 

In order to apply the iteration reconstruction algo- 
rithm described in Eq. (pS)) we need to know the over- 
all detection efficiencies T5775 and Tjrjj. In principle, 
they can be determined by knowing parameters of the 
experimental setup. However, fragility of the experimen- 
tal alignment enforces the determination from the ob- 
tained experimental data. The values of detection effi- 
ciencies can either be derived from the measured covari- 
ance between the signal and idler photoelectron num- 
bers Cs and c/ or alternatively by applying a method 
described in Sec. VI below that relies on finding the 
best fit to the experimental data. This method applied 
to the set of data (b) has provided Tsrjs = 0.207 and 
TiTji = 0.205 that have been used in the reconstruction. 
These values take into account quantum detection effi- 
ciency of the iCCD camera as well as losses occurring 
in the setup (frequency filters, reflection on the output 
plane of the crystal and mirror). The reconstruction al- 
gorithm also needs the level of dark noise that has been 
monitored in the third region-of-interest of the photo- 
cathode; Ds = Di = 0.03 (a), 0.09 (b), and 0.46 (c). 
The application of the formula in Eq. ((28)) has then 
resulted in the joint signal-idler photon-number distri- 
butions Prccins,ni) = p''°°\ns,ni) appropriate for the 
output plane of the crystal and shown in Fig. S) The 
initial joint signal-idler photon-number distribution p'"'* 
has been taken as uniform in all three cases. Covariances 
of the reconstructed joint signal-idler photon-number dis- 
tributions proc equal 90.0 % (a), 90.1 % (b), and 89.7 % 
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FIG. 3: Topological graphs of the measured histograms / FIG. 4: Topological graphs of the reconstructed joint signal- 

of signal (cs) and idler (c/) photoelectron numbers for three idler photon-number distributions Prcc{ns ,ni) for data sets 

data sets denoted as (a), (b), and (c). A black curve encircles (a), (b), and (c). A black curve encircles the area in which 

the area in which the classical inequality (I32|) is violated. the classical inequality (|32p is violated. 



(c). These numbers show the ability of the reconstruction 
algorithm to recover paired correlations that have been 
weakened during the propagation and detection process. 
The fact that the obtained covariances do not approach 
100 % has two reasons: 1) imperfect description of all 
noises occurring in the experiment, 2) numerical imple- 
mentation of the iteration reconstruction algorithm loses 
its precision with the increasing number of steps [29| . 
Both values of the parameter S given in Eq. and 
covariance C defined in Eq. ((30|) can be used for moni- 
toring convergence of the iteration process. Covariance 
C has been found to be more sensitive. Typically sev- 
eral hundreds of iteration steps are needed to arrive at 



solid (asymptotic) results as documented in Fig. [5] valid 
for data set (b). In all three cases, 10 000 iteration steps 
have been applied. 

The reconstructed joint signal-idler photon-number 
distributions proc show that the emitted fields are mainly 
composed of photon pairs that are responsible for nonzero 
covariances of the signal and idler photon numbers. Such 
fields are nonclassical in the sense that they cannot 
be described by nonnegative Glauber-Sudarshan quasi- 
distributions [37|. As a consequence, there even exist 
elements Prec{ns,ni) of the joint photon-number distri- 
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FIG. 5: Covariance C of the signal and idler photon numbers 
ns and ni as it depends on the number mt of the iteration 
step for data set (b). 



button picc that violate the classical inequality [5C 



Picc{ns,ni) < 



exp[-(ns) - (n/)]. (31) 



However, nonclassical properties manifest also in quan- 
tities which determination is based on all elements of the 
joint photon- number distribution p^ec- 
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Important nonclassical characteristics of paired 
fields 



Correlations in the signal and idler photon numbers ns 
and nj lead to narrowing of the distribution p_ of the 
photon-number difference ns — nj together with broad- 
ening of the distribution p+ of the photon-number sum 
ns + ni. The distributions p+ and are defined as 



p±{n) 



E 



E 

n5— nj —0 



{ns,ni). (32) 



The symbol S stands for Kronecker's delta. 

Fluctuations described by the distribution p_ of the 
photon-number difference can even be lower than those 
corresponding to any classical field with no correlations 
(see Fig. [5^). We then speak about sub-shot- noise cor- 
relations that can be quantified by parameter R (R < 1 
for nonclassical states): 



R ^ 



{[A{ns - ni)]') 



(33) 



Considering the experimental data, we obtain R = 0.133 
(8.7 dB) (a), R = 0.111 (9.5 dB) (b), and R = 0.117 
(9.3 dB) (c). This means that all three detected fields 
are strongly nonclassical. We note that the discussed 
narrowing can also be observed for stronger fields utiliz- 
ing correlations of photocurrents from two detectors (3l] | . 

On the other hand, the distribution p+ of the photon- 
number sum is super-Poissonian, i.e. its Fano factor F 



FIG. 6: a) Distributions p_ of the photon-number difference 
ns — ni and b) distributions p+ of the photon-number sum 
ns + ni for data set (b) (solid curves). Dashed curves give 
the distributions derived from the joint signal-idler photon- 
number distribution given as a direct product of the marginal 
signal and idler distributions; they are plotted for comparison. 



is greater than 1 [F = ((An)^)/(n)]. The suppression 
of its elements giving the probabilities of odd photon 
numbers represents its most striking feature [isj . If the 
field were composed only of photon pairs, these elements 
would have been zero. However, the presence of noise 
photons conceals this feature and so we can observe it 
only for the data set (a) obtained under low illumination 
(see Fig. |61d). 



B. Determination of photon-number statistics, 
measurement of intense fields, and role of the 
intensity profile 

The type of statistics of the emitted photon pairs is an 
important characteristic [ll, |33| • According to the theory, 
if the emission occurs in one spatiotemporal mode, the 
photon- number statistics is Gaussian (thermal). How- 
ever, the emission is usually observed into more than one 
independent spatiotemporal modes [46| and, as a con- 
sequence, the statistics of photon pairs declines towards 
a Poissonian distribution. The greater the number of 
modes, the closer the actual statistics to the Poissonian 
distribution. In the experiment the situation is more 
complex because of noises superimposed on the emitted 
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paired field. The theory presented in Sec. VI below allows 
in principle to determine the number of paired modes and 
so extract the statistics of photon pairs. When the re- 
construction algorithm is applied, we can only determine 
the statistics of marginal signal and idler fields and de- 
duce the type of statistics of photon pairs from them. 
We note that a Fano factor F is commonly used to judge 
the type of statistics, or more precisely, the declination 
of statistics from the Poissonian one. 

The Fano factor is also extraordinarily important for 
the quantification of the effect of presence of more than 
one photon in the area of one macro-pixel at the pho- 
tocathode. If the probability of having more than one 
photon at one macro-pixel is non-negligible the statistics 
of photoelectron numbers [/(cs, c/)] decline from photon- 
number statistics [precins, n-i)]- The fact that one macro- 
pixel cannot resolve photon numbers leads to a system- 
atic decrease of the Fano factor of a detected field. The 
stronger the field the smaller the value of the Fano fac- 
tor. This effect is in its nature the same as a dead-time 
effect in time-resolved detection. However, this effect can 
be corrected using appropriate transfer matrices. When 
stronger fields are measured, the transfer matrices 
and should even be corrected with respect to the 
field intensity profile as suggested in Sec. IV. 

Data set (c) has been obtained in the discussed regime 
and the regions-of-interest were composed of Ns = Nj = 
6528 macro-pixels. Here, the Fano factors of the marginal 
distributions of detected photoelectrons are Fs = 0.996 
and Fj = 1.008. We can see in Fig. [7] how the Fano fac- 
tors Fs and Fj of the marginal distributions derived from 
the reconstructed joint signal-idler photon-number distri- 
bution prcc depend on the form of transformation matri- 
ces K^''^^ and K^''^', in more detail on parameter AI 
(M = Ms = Ml) giving the number of areas inside the 
detection region. In the k-ih area there are pixels illumi- 
nated by intensities greater than (fc— 1) A/ and lower than 
/cA/, fc = 1, . . . , M [A/ = /max/M, /max being the maxi- 
mum intensity found in the profile] . The Fano factors Fs 
and Fi plotted for M = in Fig.[7]are obtained assuming 
the transfer matrices given in Eq. (ITOl) [Ns,Ni — > oo]. 
We can see from the curves in Fig.[7]that the more precise 
the form of transfer matrices the better the elimination 
of the effect and so the greater the values of Fano factors 
Fs and Fj. We can also deduce that it is sufficient to 
divide the detection region-of-interest into several areas 
to arrive at reliable results. The reconstruction of fields 
described by data set (c) clearly demonstrates the ability 
of the method to cope with this problem. 

The Fano factors of marginal distributions have been 
determined as Fs = 1.32 (a), 1.126 (b), and 1.106 (c) 
and Fi = 1.33 (a), 1.126 (b), and 1.165 (c). These values 
show that the observed down-conversion has been emit- 
ted into several tens or even hundreds of independent 
spatiotemporal modes. 

Provided that the generated joint signal-idler field is 
multi-mode and the photons are preferably generated in 
the spontaneous process, we can assume that relative 
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FIG. 7: Fano factors Fs (solid curve) and Fi (solid curve 
with *) of the marginal signal and idler photon-number dis- 
tributions obtained from the reconstructed joint signal-idler 
photon-number distribution prcc as they depend on the num- 
ber M of areas defined inside the detection region-of-interest. 

phases of different modes have random values. In this 
case, the joint signal-idler field can fully be character- 
ized by a joint signal-idler quasi-distribution of signal and 
idler integrated intensities Pw- The quasi-distribution 
Pw can be uniquely derived from the joint signal-idler 
photon-number distribution Prec using the decomposition 
into Laguerre polynomials (l2l . [37| . Due to pairwise char- 
acter of the emitted fields, the joint signal-idler quasi- 
distributions Py/ of integrated intensities attain negative 
values in certain regions [5l|, [s^l ■ 

VI. FIT OF THE EXPERIMENTAL DATA 
USING THE MODEL OF SIGNAL AND NOISE 

An alternative approach for obtaining a joint signal- 
idler photon-number distribution in the output plane of 
the crystal can be developed assuming a certain form of 
this distribution. We can assume for physical reasons 
that the overall field can be described by a certain form 
of superposition of signal and noise and can also be de- 
composed into three independent contributions. The first 
contribution describes photon pairs that are inside 
independent modes with mean photon-pair numbers hp. 
The second (third) contribution takes into account the 
presence of noise in the signal (idler) field and is com- 
posed of ms (jni) independent modes with mean photon 
numbers hs (6/). 

On the experimental side of the problem there are five 
first- and second-order moments of the measured pho- 
toelectron numbers: (eg), (c/), (cg), (cj), and (csc/). 
Moreover a reliable determination of overall detection ef- 
ficiencies Tsr\s E^nd T/jy/ in the experimental setup is dif- 
ficult. That is why we can consider the efficiencies T5775 
and T/77/ as parameters that should be determined from 
the experimental data. We thus have eight independent 
parameters and five measured quantities. The require- 
ment of minimum entropy of the joint photoelectron dis- 
tribution used for fitting the experimental data can be 
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FIG. 8: Topological graphs of a) joint signal-idler photon- 
number distribution Pfit{ns,ni) obtained by the fitting 
method and b) difference Ap{ns,ni) of photon- number dis- 
tributions pflt and prcc obtained by the fitting and recon- 
struction methods, respectively [Ap{Ns, Ni) = pmins ,ni) — 
Piccinsjni)] for data set (b). The black line identifies the 
elements pmins jni) violating the classical inequality (|32p . 



applied to find the most appropriate form of the joint 
signal-idler photon-number distribution p (details can be 
found in a forthcoming publication). 

The application of the method to data set (b) has 
provided the values Tsrjs = 0.207 and Titji = 0.205 
that were used in the reconstruction in Sec. V. Val- 
ues of the remaining parameters have been determined 
as TTip = 628, bp = 0.066, ms = 0.46, bs = 0.173, 
m/ = 0.018, and bi = 2.32. Comparison of the obtained 
joint signal-idler photon- number distribution pfn with the 
distribution prec revealed by the reconstruction is given in 
Fig. [51 We can see that the distribution pat is 'narrower', 
i.e. it contains less noise {Fs = 1.066, Fj = 1.068). 
This is in agreement with the expectation that the fit- 
ting method is by its nature more efficient in eliminating 
the noise. We note for comparison that using the fit- 
ting method covariance C of the signal and idler photon 
numbers ns and nj equals 99.7 % and R — 0.028. How- 
ever, compared to the reconstruction method, the fitting 
method does not take into account the number of macro- 
pixels and also cannot be applied to more intense fields. 



VII. CONCLUSIONS 

We have developed a method for the reconstruction 
of a joint signal-idler photon-number distribution using 
the measured histograms of photoelectron numbers and 
an iteration expectation maximization algorithm. In the 
framework of a general detection theory we have found 
formulas for the transfer matrices that give linear rela- 
tions between elements of a photon-number distribution 
and the corresponding photoelectron distribution. These 
formulas take into account finite quantum detection effi- 
ciencies, the level of dark counts as well as finite numbers 
of detection macro-pixels. Special formulas appropriate 
for very weak as well as high illumination intensities have 
been found. A method for the inclusion of a transverse in- 
tensity profile into the form of transfer matrices has been 
suggested. Three joint signal-idler photon-number distri- 
butions differing in mean photon-numbers have been re- 
constructed using the developed method. Some of their 
elements violate a classical inequality. Fluctuations of 
the difference of signal and idler photon numbers are 
highly suppressed due to pairing of photons in all three 
cases (sub-shot-noise correlations). Moreover there oc- 
curs a partial suppression of elements corresponding to 
odd photon numbers in the distribution of the sum of sig- 
nal and idler photon numbers for the weakest measured 
field. The developed reconstruction method has been 
compared to a method that provides the best fit of the 
experimental data assuming a joint signal-idler photon- 
number distribution in the form of superposition of signal 
and noise. The power of the reconstruction method to 
eliminate noise has been found weaker on one side. On 
the other side, it allows more realistic description of the 
detection process. This is invaluable for higher detector 
illumination intensities. We believe that the developed 
reconstruction method will stimulate a broader use of 
iCCD cameras as photon-number-resolving detectors. 
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Appendix A: Determination of an effective detection 
quantum efficiency 

We assume that a Poissonian field with mean photon 
number ji and statistical operator g, 

e= ^^exp(-^)|n)(n|, (Al) 

n=0 
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impinges on a detector with quantum efficiency tj and 
dark-count rate d. The probabihty p^°^^ of registering a 
photon is given as 

= l-(l-d)cxp(-77/i); (A2) 

the detection operator D has been introduced in Eq. ([3|) . 

If there is just one photon in the Fock state {g — |1)(1|) 
in a detected field, the probabihty p^°'^^ of its counting 



equals 

/°^'^ = l-(l-d)(l-r7). (A3) 

The requirement of equal detection probabilities p^""^^ 
and p^°^^ results in an effective quantum efficiency ri°^ 
depending on /i: 

77°ff(M) -l-cxp(-77M). (A4) 
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